#########################
### Results for Justice as Checks and Balances
### Geographic analysis for appendix 
##########################

### This scrips creates maps for the supplementary material
##  This code creates:
## Maps for Online Appendix (Figure A5, A6) 
## Tables of geographic variation Table A4 and A5 


# Preliminaries --
rm(list=ls())
gc()


###########################
# Relevant libraries
# Next, we download packages that H2O depends on.
pkgs <- c("plm", "dplyr", "magrittr", "ggplot2", "lmtest", "stargazer",
          "aod", "multiwayvcov",  "foreign", "sf",
          "RColorBrewer", "st"
)
for (pkg in pkgs) {
  if (! (pkg %in% rownames(installed.packages()))) { install.packages(pkg) }
  lapply(pkg, library, character.only=T )
}

#### Read agn_indios
# read agn_indios.csv
setwd("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/")
#agn_indios <- read.csv('agn_indios3.csv')
agn_indios <- read.csv("C://Users/franc/Dropbox/PhD/Edgar_Tlaxcala/AGN/agn_indios_2021.csv")


#######################
##  
##  MAPS
##
#######################

##### Load map 
gerhard <- st_read("C:/Users/franc/Dropbox/Aztec/Gerhard/Mex_terr_2016.shp")
### Create regions
gerhard <- gerhard %>% mutate(region=Region,
                              region=ifelse(Region=="Neuva Espana" |Region=="Nueva Espana", "Nueva Espana", region ),
                              region=ifelse(Region=="Nueva Galacia", "Nueva Galicia", region ),
                              region=ifelse(Region=="Neuvo Leon", "Nuevo Leon", region )
)


### Create a pallete
my_pal <- brewer.pal(7, "YlOrRd")
my_pal[1] <- colors()[245]
my_pal2 <- rep(my_pal[1],7)

my_pal3 <-  c(colors()[245], "white", "red")
# The palette with grey:
cbPalette <- colors()[1:13]

my_pal <- brewer.pal(13, "Paired")

#### Map of regions --------------------------------------------------
### Figure S5.1
regions <- ggplot() + geom_sf(data = gerhard, aes(fill=region), lwd = 0) +scale_fill_brewer(palette = "Set3")+
  coord_sf() +theme_bw() + theme(legend.title = element_blank(), axis.text.x = element_text(size=6),
                                 axis.text.y = element_text(size=6)
  )
regions
### Save 
#setwd("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/online_appendix/")
#ggsave("regions.png", regions, width=6, height=6)



#### Claim Map  --------------------------------------------------
#### Figure S5.2
#### Total claims
### Group 
claims_prov <- agn_indios %>% group_by(name, NewID) %>% summarise(claims=n())

gerhard <- gerhard %>% left_join(claims_prov)

gerhard <- gerhard %>% mutate(claims2=ifelse(claims<5, "a) Less than 5", NA),
                              order=ifelse(claims<5, 1, NA),
                              claims2=ifelse(claims>=5 & claims<50, "b) 5-50", claims2),
                              order=ifelse(claims>=5 & claims<50, 2, order),
                              claims2=ifelse(claims>=51 & claims<100, "c) 51-100", claims2),
                              order=ifelse(claims>=51 & claims<100, 3, order),
                              claims2=ifelse(claims>=100 & claims<150, "d) 100-150", claims2),
                              order=ifelse(claims>=100 & claims<150, 4, order),
                              claims2=ifelse(claims>150 , "e) More than 150", claims2),
                              order=ifelse(claims>150 , 5, order)
)

####
claims_map <- ggplot() + geom_sf(data = gerhard[order(gerhard$order,decreasing=T),], aes(fill=claims2), lwd = 0) +
  scale_fill_brewer(palette = "YlOrRd")+
  coord_sf() +theme_bw() + theme(legend.title = element_blank(), axis.text.x = element_text(size=6),
                                 axis.text.y = element_text(size=6)
  )
claims_map
#### Save
# setwd("C:/Users/franc/Dropbox/PhD/dissertation/paper_colonialism_liberalims/tex/2019_version/online_appendix/")
# ggsave("claims_tot.png", claims_map , width=6, height=6)


#######################
##  
##   Tables by region
##
#######################
### Tables S5.1 and S5.2 
### Count of claims by region
table(agn_indios$region)

## For the rest of the analysis we need to delete the provinces of Valladolid and Guachinango
## Because the location is ambigous (Two Valladolids and two Guachinangos)
#agn_indios$X[duplicated(agn_indios$X)]
agn_indios_sub <- agn_indios %>% filter(name!="Valladolid", name!="Guachinango", name!="Guanchinango")
enc_mean <- mean(agn_indios_sub$per_encomienda, na.rm=T) # 0.1529982


#### NE ##################################
model_NE <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
               data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                       agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE"
               ),])
model_NE_se <- cluster.vcov(model_NE, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE"),]$name)
model_NE_se  <- coeftest(model_NE , model_NE_se)

#### Pop dec
model_NE_dec <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                   data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                           agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" & agn_indios_sub$pop_ratio<1
                   ),])
model_NE_dec_se <- cluster.vcov(model_NE_dec, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                   agn_indios_sub$region=="NE"& agn_indios_sub$pop_ratio<1),]$name)
model_NE_dec_se  <- coeftest(model_NE_dec , model_NE_dec_se)

#### Pop inc
model_NE_inc <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                   data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                           agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" & agn_indios_sub$pop_ratio>=1
                   ),])
model_NE_inc_se <- cluster.vcov(model_NE_inc, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                   agn_indios_sub$region=="NE"& agn_indios_sub$pop_ratio>=1),]$name)
model_NE_inc_se  <- coeftest(model_NE_inc , model_NE_inc_se)


##### LE strong
model_NE_str <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                   data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                           agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" & agn_indios_sub$per_encomienda>=enc_mean
                   ),])
model_NE_str_se <- cluster.vcov(model_NE_str, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                   agn_indios_sub$region=="NE"&  agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_NE_str_se  <- coeftest(model_NE_str , model_NE_str_se)

#### LE Weak
model_NE_wk <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                  data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                          agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" & agn_indios_sub$per_encomienda<enc_mean
                  ),])
model_NE_wk_se <- cluster.vcov(model_NE_wk, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                 agn_indios_sub$region=="NE"&  agn_indios_sub$per_encomienda<enc_mean),]$name)
model_NE_wk_se  <- coeftest(model_NE_wk , model_NE_wk_se)


##### Conditioned on decrease
##### LE strong
model_NE_dec_str <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                       data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                               agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" & agn_indios_sub$pop_ratio<1 & 
                                                 agn_indios_sub$per_encomienda>=enc_mean
                       ),])
model_NE_dec_str_se <- cluster.vcov(model_NE_dec_str, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                           agn_indios_sub$region=="NE"& agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean),]$name)
model_NE_dec_str_se  <- coeftest(model_NE_dec_str , model_NE_dec_str_se)

#### LE Weak
model_NE_dec_wk <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                      data=agn_indios_sub[which(agn_indios_sub$year>1592 &
                                              agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE" &
                                                agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean
                      ),])
model_NE_dec_wk_se <- cluster.vcov(model_NE_dec_wk, ~ agn_indios_sub[which(agn_indios_sub$year>1592 & agn_indios_sub$name!="Mexico" & 
                                                                         agn_indios_sub$region=="NE"&  
                                                                           agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean),]$name)
model_NE_dec_wk_se  <- coeftest(model_NE_dec_wk , model_NE_dec_wk_se)


### Units 
##### Number of Units by regression
all   <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                   agn_indios_sub$name!="Mexico" & agn_indios_sub$region=="NE")]))

dec  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                  agn_indios_sub$name!="Mexico" & agn_indios_sub$pop_ratio<1 &agn_indios_sub$region=="NE")]))

inc  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                  agn_indios_sub$name!="Mexico" & agn_indios_sub$pop_ratio>=1 &agn_indios_sub$region=="NE")]))

strg  <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                   agn_indios_sub$name!="Mexico" & agn_indios_sub$per_encomienda>=enc_mean &agn_indios_sub$region=="NE" )]))

weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                  agn_indios_sub$name!="Mexico" & agn_indios_sub$per_encomienda<enc_mean &agn_indios_sub$region=="NE" )]))

dec_str <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                     agn_indios_sub$name!="Mexico" & agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda>=enc_mean &agn_indios_sub$region=="NE")]))

dec_weak <- length(unique(agn_indios_sub$name[which(agn_indios_sub$year>1592 &
                                                      agn_indios_sub$name!="Mexico" & agn_indios_sub$pop_ratio<1 & agn_indios_sub$per_encomienda<enc_mean &agn_indios_sub$region=="NE")]))

unit_line <- c("Units", all, dec, inc, strg, weak, dec_str, dec_weak)


#### Table S5.1
#### With Robust errors
stargazer(model_NE ,model_NE_dec, model_NE_inc, model_NE_str, model_NE_wk, model_NE_dec_str, model_NE_dec_wk,
          se   = list(model_NE_se[,2] ,model_NE_dec_se[,2], model_NE_inc_se[,2], 
                      model_NE_str_se[,2], model_NE_wk_se[,2], model_NE_dec_str_se[,2], model_NE_dec_wk_se[,2]
          ) ,
          p    = list(model_NE_se[,4] ,model_NE_dec_se[,4], model_NE_inc_se[,4], 
                      model_NE_str_se[,4], model_NE_wk_se[,4], model_NE_dec_str_se[,4], model_NE_dec_wk_se[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("All", "Population", "Local Elites", "Population Decrease"
          ),
          column.separate = c(1,2,2,2),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y","Y"),
                         c("Time FE", "Y", "Y", "Y", "Y","Y", "Y","Y")
          ),
          label="dictionary_NE",
          title =("Favorable Ruling, Topics, and Balance of Powers (New Spain)") 
          #out="model_NE.tex"
          )




##### North fromtier NG and NV  ------------------------------------------------------------------
#### Table S5.2
per_enc_north <- mean(agn_indios$per_encomienda.x[(agn_indios$region=="NG"|  agn_indios$region=="NV")], na.rm=T)

#### NE ##################################
model_North <- lm(protect_new2 ~tierras_d+ abusos_d +impuestos_d + comunidad_d+as.factor(decade)+name,
                  data=agn_indios[which(agn_indios$year>1592 &(agn_indios$region=="NG"|  agn_indios$region=="NV")
                  ),])
model_North_se <- cluster.vcov(model_North, ~ agn_indios[which(agn_indios$year>1592 &(agn_indios$region=="NG"|  agn_indios$region=="NV")),]$name)
model_North_se  <- coeftest(model_North , model_North_se)


unit_north <- length(unique(agn_indios[which(agn_indios$year>1592 
                                                 &(agn_indios$region=="NG"| agn_indios$region=="NV")),]))

##### Stargazer North frontier 
stargazer(model_North , 
          se   = list(model_North_se[,2]
          ) ,
          p    = list(model_North_se[,4]
          ),
          dep.var.labels= c("Favorable Ruling"),
          covariate.labels =c("Land Conflicts","Mistreatment",  "Taxes", 
                              "Community"
          ),
          column.labels = c("All", "Population", "Local Elites", "Population Decrease"
          ),
          column.separate = c(1,2,2,2),
          object.names = F,
          #model.names = c("", "Decrease", "Increase", "Strong", "Weak", "Strong", "Weak"),
          keep = c("tierras_d", "abusos_d", "impuestos_d" , "comunidad_d"), 
          keep.stat = c("n", "adj.rsq"),
          add.lines=list(c("Province FE", "Y", "Y", "Y", "Y", "Y", "Y","Y"),
                         c("Time FE", "Y", "Y", "Y", "Y","Y", "Y","Y")
          ),
          label="dictionary_North",
          title =("Favorable Ruling, Topics, and Balance of Powers (Northern Frontier)") 
#          out="model_NE.tex"
          )

